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Abstract 



Semi-elliptic stochastic differential equations (SDEs) are common mod- 
els among practitioners. However, the value functions and sensitivities 
of such models are described by degenerate parabolic partial differential 
equations (PDEs) where the existence of regular global solutions is not 
trivial, and where densities do not exist in spaces of measurable functions 
but only in a distributional sense in general. In this paper, we show that 
Ci . for a related class of such equations regular global solutions can be con- 

structed. Moreover, the solution scheme has a probabilistic interpretation 
where the existence of regular densities on certain subspaces of the state 
space can be exploited. Prominent examples of models of practical inter- 
est belonging to this class include factor reduced LIBOR market models 
and Cheyette models. Moreover, factor reduced SDEs originating from a 
full factor model are in the class to which our theorem applies. 
^— ^ ] The result is also of interest for the theory of degenerate parabolic 

■ equations. A more detailed analysis of numerical and computational is- 

sues, as well as quantitative experiments will be found in the second part. 
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1 Introduction 



In this paper, we show that for a class of semi-elliptic equations with relatively 
rough initial data on a linear subspace (where a hypo-cllipticity condition holds), 
and regular data on the complementary linear subspace, regular global solutions 
can be constructed. Moreover, we derive a solution scheme which has a proba- 
bilistic interpretation such that it can be used to build a simulation scheme of 
the associated SDE. The class of equations considered is particularly relevant 
for valuation models derived in finance. 

1.1 Main Theorem and its Assumptions 

The prototypal application of the main theorem of this paper is the valuation of a 
financial derivative: put simply, the valuation of a financial derivative consists of 
a model of some undcrlyings (usually given by a stochastic differential equation 
(SDE)) and a value function representing the value of the financial derivative as 
a function of the modeled underlyings. In its simplest case, the value function is 
known at some future time as a function of the underlyings, the payoff function, 
and the value function of the financial derivative at earlier times can then be 
represented as some conditional expectation of the payoff function. 1 

Our main theorem makes an assumption about the relation of the model 
SDE and the payoff function, namely 

• the model SDE is semi-elliptic, i.e. may have a degenerated diffusion 
matrix, and 

• the payoff function may be rough ( i.e. in L p ) in some state variables, and 

• the subspace of smoothness of the payoff function fits a subspace where 
the semi-elliptic operator corresponding to the SDE is not smoothing, i.e. 
is not even micro- hypoelliptic. 

In Section 2.1.1 we illustrate that this assumption is a natural situation in 
derivative pricing models. 

1.2 Layout of the Paper, Proof of the Main Theorem 

Our main theorem is motivated by stochastic differential equations and payoff 
functions encountered in the realm of valuing financial derivatives. We con- 
sider semi-elliptic models which satisfy a hypo-cllipticity condition on a linear 
subspace and have adapted payoffs with mixed regularity. 

Our main theorem is an extension of Hormander's well-known result concern- 
ing a sufficient condition for the hypoellipticity of linear differential operators 
with variable coefficients (the so called Hormander condition, cf. [12]). Smooth- 
ness (i.e. C°°-regularity) of the density of the associated SDEs is a subject of 
Malliavin calculus. Indeed a Malliavin-type estimate obtained in [23] is crucial 
for our main theorem. 

Moreover, our proof consists of a constructive PDE-scheme. This scheme 
can be turned into a Monte-Carlo scheme for the original associated SDE. 

1 The more general case where intermediate payoffs depend on future values of the derivative 
is usually a straight forward extension, defining a backward-induction of value functions. 
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The paper is structured as follows: Section 2 gives a thorough motivation 
for our considered class of SDEs and payoff functions and discusses their re- 
lationships to PDEs, transition densities and the calculation of sensitivities in 
the context of semi-elliptic models. Section 3 gives a list of some models where 
our main theorem readily applies. These models are common among practition- 
ers. Section 4 extends our considerations to general reduced diffusion models. 
Section 5 reviews the existence of global weak solutions from the stochastic per- 
spective and motivates the main theorem from a mathematical perspective. In 
Section 6 we state and prove our main theorem formulated as a global exis- 
tence result of the associated PDE. Section 8 briefly sketches the application 
of our constructive proof from Section 6 to the Monte-Carlo simulation of the 
corresponding SDE. Section 9 concludes the paper with some further remarks. 

2 Motivation 

2.1 Semi-Elliptic Stochastic Differential Equations in Fi- 
nance 

Semi-elliptic stochastic differential equations (SDEs) are common models among 
practitioners. One approach to obtaining such a model is to start from a model 
SDE with full rank diffusion matrix (i.e., full factor model) and perform a princi- 
pal component analysis (PCA) to build a factor reduced model with a low-rank 
diffusion matrix. The result is a 'reduced model' which uses only a reduced set 
of, e.g., Brownian drivers while retaining the full state space dimension. The 
convection term, i.e. drift, is not affected by the reduction. A prominent exam- 
ple of a model where a factor reduction is usually applied is the LIBOR market 
model. A factor reduction can be advantageous in terms of calculation speed, 
see [14] and of course noise reduction. Often the model is reduced to factors 
(i.e., modes) which are of particular importance to the application and PCA 
is just one way of determining important factors. If, for example, a model is 
used as a pricing model, one might determine important factors from product 
properties. Practitioners usually attribute importance to factors like level, slope 
and curvature, see [19]. 

Another approach to obtaining a factor-reduced SDE is to use a specific 
modeling assumption leading to a low factor model with, e.g., low Markov di- 
mension. A prominent example here is the Cheyette model: starting from an 
infinite dimensional HJM model a specific (separable) instantaneous volatility 
structure will lead to an, e.g., two dimensional Markovian one factor model. 
These examples belong to a general class of stochastic processes X t with values 
in W 1 and satisfying a stochastic differential equation 

dX(t) = b(t,X(t))dt + a(t,X(t))dW(t), < t < T, X(0) = x. (1) 

Elementary SDE theory tells us that the global existence of Levy-continuous 
solutions can be derived if the coefficient functions satisfy a global Lipschitz 
condition (cf. discussion in Section (5) below). The Feynman-Kac formalism 
shows that for a considerable class of data functions / (related to the payoffs in 
finance) and stochastic processes X the function 

(t,x)->E(f(X t )) (2) 
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and its derivatives (related to derivative prices and Greeks in finance) are de- 
termined by the solution of a Cauchy problem on [0,T] x W l (where T > is 
an arbitrary finite horizon) 

(3) 

u(0,x) = f{x), 
where 

«=E»%, (4) 

j=l 3 

are vector fields with < i < m. Hormandcr provided a sufficient condition for 
the existence of regular densities. However, in the case of (1) this Hormandcr 
condition is satisfied only on a certain subspace of M. n in general. In most 
applications in finance, reduced models are of a special form where the vector 
fields (4) live in a fixed subspace such that the Cauchy problem (3) has a certain 
block structure (we list a set of the most popular models below) . For this reason 
we emphasize this case because it is of particular interest for the practitioner. 
It turns out that our main result can be extended to a more general result 
which does not depend on the block structure and which is an extension of 
Hormander's result. We state this generalization (the proof method is similar). 

2.1.1 Relation of the Payoff and the Model 

A factor reduction, i.e. the reduction of the brownian drivers, is often motivated 
by computationally feasibility: the model should be computational feasible and 
still provide an accurate representation of the perceived underlying dynamic 
(e.g., historically observed "factors"). However, in addition the modeling has to 
consider the application. For the valuation of a payoff, the degree of regularity 
of the payoff function will give weight to certain aspects of the underlying's 
dynamics. This weighting of "factors" is often different from the distribution of 
historically observed factor dynamics. It is exactly this weighting which leads 
to the fact that a meaningful semi-elliptic model for financial derivative pricing 
will most likely fall into the class of models considered by our main theorem. 2 

Example A simple example of the situation described is the modeling of an 
interest rate curve and the valuation of an (interest rate) spread option. The 
interest rate curve can be formalized as an infinite dimensional space. One 
approach to modeling it is to focus on a discrete set of interest rates (e.g., 
forward rates with fixed period length) and define the interest rate curve using 
appropriate interpolations. Let us consider the model 

dL k (t) = f i k (t)dt + a(t)dW 1 (t)+ 7 (t)dW 2 (t), 
dL k+1 (t) = fi k+1 (t)dt + a^dWxit) - 7 {t)dW 2 {t) 

for two adjacent forward rates of fixed period length L k and Lk+i- Assume that 
W\ and W 2 are independent Brownian drivers. In this model a is the instanta- 
neous volatility of a joint interest rate movement while 7 is the instantaneous 

2 We will give a selection of some common models in Section 3. 



5 



volatility of the spread. In other words: a is the loading factor of changes in the 
curve level and 7 is the loading factor of changes in the slope. From empirical 
data it appears to be reasonable that 7 is much smaller than a. Hence a mod- 
eler may favor a computationally more efficient approximation of the model by 
neglecting 7 and simply consider the reduced model 

dl k (i) = fx k (t,L k ,L k+1 )dt + a(t)dW(t), 
dL k+1 (t) = fx k+1 (t,L k ,L k+1 )dt + a(t)dW(t) 

Let T = + AT denote a small time step. Choosing the time T-terminal 
measure, i.e., the numeraire N(T) = 1 and the associated equivalent martingale 
measure Q we have: for a payoff function / = f(L k (T), L k+ \(T)) the value 
function is given by the expectation operator E Q (f(L k (T), L k +i(T)) \!Fo)- 

Consider the valuation problem of a spread option with payoff 

h{L k {T),L k+1 {T)) := max(L fc (T)-L fc+1 (T),0). 

While the historical movement of the interest rate curve is such that a larger 
part of its dynamic is captured by a one-factor model, the reduced model is not 
suitable for valuing the spread option. The reduced model completely ignores 
the spread volatility. It is 

d(L k (t) -L fc +i(t)) = {^ k+x {t)L k+1 {t)-ii k {t)L k {t))dt + 2 1 dW{t). 

For the value function we find 

E{h{L k {T),L k+1 {T)) |Jo) - max(L fe (0)-T fc+1 (0),0)+7O(VT), 

as T becomes small. Compare this to a linear payout function of the spread, 
say 

f 2 (L k (T),L k+1 (T)) := c(L k (T) - L k+l {T)) 

then we indeed may neglect 7 since ryO(VT) does not enter the value function. 
We have (as T becomes small) 

E(f 2 (L k (T),L k+1 (T)) |To) = c(L k (0)-L k+1 (0)) + O(T). 

In other words, there is a relation between modeling and payoff. The relation can 
be described as follows: A model has to consider a source of risk (i.e. diffusion) 
in states where the payoff is discontinuous or strongly convex, and, the model 
may neglect a source of risk (i.e. diffusion) where the payoff is smooth and 
only weakly non-linear. This aspect of modeling corresponds to the regularity 
assumptions formulated in our main theorem. 



2.2 Relation to PDEs, Existence of Densities 

While a factor reduction is straight-forward on the SDE level and easy to imple- 
ment, it may severely affect the sense in which the solution of the corresponding 
PDE, i.e., the density exists. The most important situation we consider in this 
paper (at least from the point of view of the practioner) is the following: Con- 
sider a matrix function x — > (a,ji) d,m (x), 1 < j < d on 1" with n > d, and m 
smooth vector fields of dimension d 

d d 



G 



where < i < m. Consider an additional vector field of dimension n — d 

n „ 

B := V „ ,6) 

j=d+l J 

Consider the Cauchy problem on [0, T] x R n (where T > is an arbitrary finite 
horizon) 

^ = \YZi^+{A + B)u 

(7) 

u(0,x) = f(x). 

Remark 1. The Cauchy problem (7) is written in the block structure form 
which we have in most cases of practical interest. Generalisations which are 
coordinate-independent are considered in section 4 and in section 7. 

Remark 2. Note that the Cauchy problem (7) may be written in the more 
familiar form 

du 1 sr^k n * d 2 u sr^n Q u ,a 

dt 2 l^i. j I "i.l ■■• l^i=l Mi dx* ~ U ' 

(8) 

u(0,x) = f(x), 

for some k < n, and where (a*j) is a diffusion matrix and [ii are drift terms 
defined in terms of coefficients of the vector fields Ai. Note that we have k < d. 
If (a*j ) is uniformly elliptic on the subspace M. k , then we are back in the situation 
of [9]. 

Assume that (7) satisfies the Hormander condition with respect to the sub- 
space M. d , i.e. assume that 

{A h [Aj,Ak] , [[Aj,A k ] , Ai] , ■ ■ ■ |1 < i < m, 0<j,k,l---< m} (9) 

spans M. d at each point x. Here [., .] denotes the Lie bracket of vector fields. 
If n = d we are back in the situation of Hormander's class of hypoellipticty. 
The existence of regular solutions of the Cauchy problem in this case is well 
known (cf. [12] ). Indeed Hormander's result shows us that a family of smooth 
transition densities exists if (9) holds for every x € W 1 (where n = d). This is 
no longer true if n > d. In this paper we show that in this case we still have 
regular solutions if the Hormander condition holds on a subspace where data 
are locally measurable (i.e. L^ c ) and here data are smooth (and satisfying some 
growth condition) on the complementary space where the Hormander condition 
does not hold. 

Note that a consequence of a factor reduction is that the density may lose 
regularity. Often a classical (regular) density will no longer exist. Consider the 
simple example in the following section 2.3. 

2.3 Example: A simple two dimensional toy model 

Let us consider the following Cauchy problem on [0,T] x K 2 : 

du 1 „2 d 2 u _j_ , , du 

dt — 2 a ~r P dx 2 ' 

(10) 

u(0,x) = f(xi) +g{x 2 )- 
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The solution of this equation is 



1 / (xi—yi) 



'2-Kt 



f(yi)— 



This leads us to a 'distributional density' of the form 

p(t,x,y) := exp (- ^"^ ) S(x 2 + fit - y 2 ). (12) 

Indeed, formally we have (let us denote dy = dy\dy 2 ) 
I (f(yi)+9(V2))p(t,x,y)dy = 

J K /(2/i)^=t exp (- ^"a^ j <5(z2 + M* ~ V2)dy 1 dy 2 
+ J R g(^2)^|s- exp (- (xi 2 ~3t } ) <K X 2 + /xt - y2)dyidy 2 ^ 



= /b /(yi) exp (" '^Jt^ ) + 9{x 2 +fd- y 2 ) 
= u(t,x), 

assuming that some 'kind of Fubini theorem' can be applied. We see from this 
example that no regular density exists except in a distributional sense. This 
simple example shows that no regular density exists in a situation of degenerate 
diffusion models. More particularly, the regularity theory of densities in the 
context of Malliavin calculus does not apply. Note that the associated process 
is X = (X!,X 2 ) with 

dXi = adW 

(14) 

dX 2 = [idt. 

Such a process exists as is well-known in the context of stochastic analysis 
(cf. citation of theorem 1 below). The simple example shows that no regular 
density exists in a situation of degenerate diffusion models with n > d. On the 
other hand, (14) (and thus (10)) could just be the result of a factor reduction 
of the process 

dXi = adW 

(15) 

dX 2 = a 2 dW 2 +ndt. K ' 

where the factor reduction (i.e., dropping dW 2 ) may have been applied because 
c 2 was considered to be small. 

2.4 Discussion of the Cauchy-problem (7) from the per- 
spective of Malliavin calculus, microhypoellipticity, 
and PDE-methods 

The example above shows that operators of the form (7) are not hypoclliptic 
for B ^ in general. Recall that a differential operator L with C^-coefficients 
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is called hypoelliptic on an open set ft C R™ if for any distribution u on Q, is in 
C°° if Lu^C°°. Now, for u with 

u(t,x) = [ /(2/1)— =^ cxp (- <yXl ^ \ dy x + g(x 2 + fit). (16) 



we have 

<9u 1 2 <9 2 u <9u 

where 3 may be C , and this shows that the operator in example 1 is not 
hypoelliptic. This provides us with a rigorous argument that there is no regular 
density for the Cauchy problem in (10), because it is known that for linear 
operators with constant coefficients hypoellipticity is equivalent to the existence 
of a fundamental solution in C°° on R™ \ {0}. 

It is of interest here to note that operators with constant coefficients may not 
be hypoelliptic even if they are obtained by freezing coefficients of hypoelliptic 
operators with variable coefficients. For example the hypoelliptic operator on 
R 3 of the form 

d 2 ( 9 d\ 2 



dx 2 \ dy dz 
has a counterpart (c = some constant) 



d 2 fd d\ 2 



9a; 2 \ dy dz 

which is not hypoelliptic, because h(cy — z) is a solution of L c h = for a C 1 
function h. Note that Hormander's result is sensitive to such differences, since it 
is concerned with microhypoellipticity which implies hypoellipticity. Consider 
equation (7) with B = 0, i.e. consider 

(20) 

u(0,x) = f(x). 

If Hormander's condition is satisfied (20), then the operator of (20) is microhy- 
poelliptic, i.e. 

WF(u) = WF(iu), (21) 

where WF(«) denotes the wave front set of a distribution u, i.e. the intersection 
of the characteristic varieties of pseudo-differential operators P of order zero 
which satisfy Pu £ C°° . Indeed, it is clear that the operators of (7) are not 
microhypoclliptic for B ^ in general. This is reflected by the fact for pseu- 
dodiffcrential operators of negative order the characteristic variety equals the 
whole cotangential bundle T°f2 of the domain ft. Hence, the microhypoclliptic 
theory is designed for pscudodiffercntial operators of order (at least), and this 
does not apply in our case. 

From the probabilistic point of view the Cauchy problem (7) has a Levy- 
continuous solution (cf. the discussion in section 2 below). For smooth data 
we immediately get regularity of the first order time derivative and a classical 
solution (cf. remark 3 below). However, especially in finance, we are interested 
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in less regular data (payoffs). The regularity theory of densities in the context 
of Malliavin calculus does not apply directly (as we remarked in [9]). However, 
it may still be useful, especially for the analysis of the subspace of dimension d 
where the Hormander condition holds. Let us consider (7) from a probabilistic 
perspective. The associated Stratonovic integral of a process starting at x E W n 
is 

pt »« ft 

/ A {X s )ds + V / A k (X s )odW k (s) (22) 
Jo k=1 Jo 

where Wi denote Brownian motions and o indicates that the integral is in the 
Stratonovic sense. If n = d and the Hormander condition holds, then the 
associated covariance matrix process at is defined by 



X t = x 



°t = z- 1 



t 



Zt 1,T (23) 



Z s A s (X s )Al {X S )Z A S 
where Z is a matrix- valued invcrtible process defined by 

i(X s )doWi(s), (24) 



pt " pt 

I d - / Z s DA (X s )ds- V / Z s DAi( 
Jo l=1 Jo 



and Id denotes the ci-dimensional identity matrix. From the perspective of 
financial applications it is interesting to note that at is almost surely invcrtible, 

a- 1 G ^ (25) 

for every real number p and t in some arbitrary finite time horizon [0, T], because 
this property is shared by many stochastic volatility diffusion models. The 
reduced models where it does not hold but holds only on a linear subspace are 
considered in this paper. 

The theory of partial differential equations has investigated systems of equa- 
tions of hyperbolic-parabolic type for a long time. One of the crucial conditions 
for such systems is the so-called Kawashima-condition. Different equivalent 
forms of this condition are considered in [24]. Several equivalent conditions of 
global existence and regularity have been developed for systems of the form 

n n 

D°w t + Dj w Xj - Ejkw ^ k + Fw, (26) 

j=l j,k=l 

where Di, E° , and F are are m x m-matrices. For general results these matrices 
are considered to be constant. In this case, for intial data 

w e H s (R n ) n V (R n ) , (27) 

and if D j , l<j< n and E' k , l<j,k< n are real symmetric matrices, and 

jk 

for £ G S 71 " 1 (the n— 1-dimcnsional sphere in W 1 ), then for dissipative systems 
there are constants 5, C > where we have 

\\D l x w(t)\\ L > < C Ml^ojl + (1 + <F +i/2 |M|J (29) 
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for all integers / < s. For variable coefficient matrices the existence results are 
weaker (additional conditions, small data or less regularity). Observe that such 
results are derived for data with uniform regularity assumptions. Our theorem 
is different and special in this respect. 

2.5 Relation to Sensitivities (Greeks) 

In our main theorem we will construct a regular solution to the associated 
(backward) PDE of a factor reduced SDE, given that 

• the regularity of the initial condition of the PDE, i.e. the regularity of 
the associated payoff function, is (in a certain sense) compatible with the 
diffusion matrix, namely that we have smoothness of the initial condition 
where we lack diffusion. 

Note also that this assumption is important, or even a requirement, when 
computing Greeks. If, for example, the payoff function is non-differcntiable 
in a coordinate for which the underlying SDE lacks diffusion the correspond- 
ing (second order) partial derivative, i.e. the sensitivity or Greek, will be not 
unique (or unbounded). In this situation, the "ignorant" numerical calculation 
of path-wise Monte-Carlo Greeks will not converge and the numerical calcu- 
lation of a Greek involving likelihood rations cannot be performed due to the 
lack of a (regular) density. The constructive scheme involves regular densities 
on the subspace where the diffusion lives and these densities can be used when 
calculating Greeks numerically. Thus, the constructive iteration scheme used in 
the proof of our main theorem may itself result in a numerical (Monte-Carlo) 
scheme to calculate robust Monte-Carlo Greeks in the setup of a reduced factor 
model. 



3 Applications from Finance 

3.1 Factor Reduced Lognormal LIBOR Market Models 

Reduced factor (LIBOR market) models are common among practitioners. In 
[10] a numerical scheme of a reduced factor LIBOR market model is given by 
considering the full proxy scheme on a linear subspace; in [7] by considering 
the partial proxy scheme. Both schemes require knowledge of transition proba- 
bilities. In [10] it was possible to give the transition density for the full factor 
model (WKB expansion) , while for the reduced factor model only the transition 
density of an approximating scheme was considered. As an important example, 
let us look at the LIBOR market model with tenor structure < T\ . . , < T n+ \ 
in terminal measure P n +i (induced by the terminal zero coupon bond B n+ i(t)). 
The dynamics of the forward LIBORs Li(t), defined in the interval [0,Tj] for 
1 < i < n, are described by 

du = " 8 3 ul gtn &t + L Jdw(d) 

4ti 1 + 6jLj (30) 
=: rfi^^U + L^J&W^, 

where 5i = T i+ i —Ti are day count fractions and t —> ji(t) = (7i,i(t), . . . , 7i,d(i)), 
< t < Ti, are deterministic, bounded and smooth volatility vector functions. 
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We denote the matrix with rows by T and assume that T is invcrtible. For 
simplicity, assume that in the following T(t) = T does not depend on t (it is 
just a technical matter to generalize to time-dependent coefficients) . In (30), 
(W^ n+1 \t) \ 0<t<T n ) is a standard d-dimensionaJ Wiener process under the 
measure P n +i with d, 1 < d < n, being the number of driving factors. First, we 
consider the full-factor LIBOR model with d = n in the time interval [0, T{). 
The standard log-transformation leads to Ki := lnL^, 1 < i < n, and 

dKl = j-dL l — ^jd(Li) = 

+ e Kl , • ■ • , e*»)) dt + -fJdW^+V (31) 

= : ^ dt + dW^K 

where fi L is given in (30). In these coordinates the associated infinitesimal 
generator is related to a strictly parabolic operator, while in the original co- 
ordinates parabolic degeneracies appear. Nevertheless, parabolic degeneracies 
appear even in logarithmic coordinates Ki for reduced factor models. In such 
models one considers d driving factors for some d, < d < n instead of n driving 
factors up to time T\. Then, instead of an invertible matrix T, we have a matrix 
F typically of rank d, e.g. the factor-reduced LIBOR market model typically 
used by practitioners takes the form 

dK = n K (t,K)dt + FdW, (32) 

where F is a n x d matrix with d < n (e.g. 40 x 5 matrix). Reduced factor 
models arc popular among practitioners, because they are computationally par- 
simonious, while the essential dynamical features of the LIBOR market model 
can be preserved if d is not too small. Note that 

FF T (33) 

is a n x n matrix of rank d while F T F is a d x d matrix of rank d. If G = 
(fd+i, ■ ■ ■ fn) is the matrix consisting of the eigenvectors fi of ker (FF T ), then 
the system (32) is equivalent to 

d (F T K) = F T fi K (t, K)dt + F T FdW( n +V 

(34) 

d(G T K) = G T fi{t,K)dt 1 

because G T F = 0. The evaluation problems for options and Greeks in this 
type of model lead us directly to degenerate parabolic Cauchy problems (cf. 
(35) below). The interest in regular global solutions to projective degenerate 
equations is not limited to reduced factor LIBOR market models, of course. So 
let us formulate it in a fashion which looks less specific, i.e. without reference to 
lognormal coordinates, and as the more familiar initial-value problem. Instead 
of specific loading factor matrix T we have a volatility matrix aa T , which may 
depend on time and assets. In logarithmic coordinates this leads to the Cauchy 
problem 

§ + \ Ei,j=i (™"% (t,*) + £?=i MM) & = o, 

(35) 

u(0,x) = f(x). 
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with < k < n to have global regular solutions. For simplicity of notation, we 
use the abbreviation 

Oij(t,x) := i (ffcr T ) y . (t,x), (36) 

where we keep in mind the fact that we assume that for each (t, x) in the domain 
of (35) the matrix aij(t,x) has a decomposition as in (36). Note that condition 
(36) is not trivial; a sufficient condition is that all a.y S C 2 are bounded and have 
bounded first and second order derivatives. This assumption is important if we 
want to make the link to probabilistic schemes and to Monte-Carlo methods, as 
we do in the second part. However, for the statement of the main theorem we 
do not need these assumptions. From a practical point of view this seems not 
to be a serious restriction. 



3.2 Low Dimensional Markovian HJM Models, Cheyette 
Model 

The LIBOR market model can be interpreted as a finite dimensional Markovian 
Heath- Jarrow-Morton (HJM) model with a (usually) high Markov dimension 
n. We may observe the situation of semi-elliptic equations in lower dimensional 
Markovian HJM models, e.g. as for the Cheyette model. The Cheyette model, 
cf . [3] , can be derived from an HJM model by assuming a specific structure of the 
instantaneous volatility structure. The volatility structures were first considered 
by Ritchkcn and Sankarasubramanian (1995), cf. [22]. For an introduction to 
the HJM framework see, e.g., [6]. 

Given an HJM dynamic for the instantaneous forward rate 

d/(i, T) = a(t, T) dt + a{t, T) dW(t), 

(37) 

/(0,T) = ./o(T) 

endowed with a special volatility structure 

a(t,T) := g(t) h(T), (38) 

where g : [0, T*] h-> K\{0} denotes a deterministic function and h : [0, T*} x Q, h-> 
M. m an m-dimensional Markov process. Then the short-rate process is given by 

r(t) = /(0, *) + *(*), (39) 

where 

dX(t) = (Y(t) - n(t)X(t)) dt + rj(t) dW(t), X(0) = 0, 
dY(t) = (?y 2 (t) - 2n(t)Y(t)) dt, Y(0) = 0, 

and 

r)(t) = a(t, t) = g(t) h(t), K(t) =-^§-- 

git) 

Here, the factor reduced model (40) results from a modeling assumption and 
not from a PCA factor reduction. 
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3.3 Heston Model 



The Heston stochastic volatility model is given by 

dS(t) = (iS{t)dt + y/v(fjS(t)dWi{t) 
dv = k(6- u(t))dt + Cv / ^)dty 2 (i) 
dWi(t) dW 2 (t) = pdt 

The numerical simulation of the Heston model is known as a challanging prob- 
lem. Exact simulation is possible, but requires computationally expensive Fourier 
transforms, see [1]. Here, a proxy simulation scheme [8] can be an alternative, 
using a simple Euler scheme for path generation and to correct the probability 
densities. 



3.4 SABR Model 

The SABR model [11] is given by 

dL(t) 
da 

dW^t) dW 2 {t) 

where L is some forward rate. 3 The 
factor LIBOR market model with sm 



a(t)L{tfdWi{t) 

a(t)a(t)dW 2 (t) 

pdt 

odel may be used to endow a reduced 
modeling. 



3.5 Stochastic Volatility CEV LIBOR Market Model / 
LMM-SABR Model 

A popular extension of the LIBOR Market Model is to endow the model with 
SABR-like stochastic volatility CEV dynamics, see [21]. One such extension is 

dL, = pf(t,L)dt + if 
dX = a\dW d+1 . 

Using the transformation 

K . . . / L]^' for ft ^ 1 
1 ■ \ ln(Li) for ft = 1 

we can write the model with a block diffusion structure as 

dA = a\dW d+1 , 
d{F T K) = F T p K {t,L)dt + XF T FdW {d \ 
d(G T K) = G T p K (t,L)dt, 

where F, G are as in Section 3.1. 

3 The model may also be used to model other underlyings. 
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4 General reduced Diffusion models 



Although the examples above show that there is a considerable class of reduced 
market models which can be put into the form (7), a more general form of a 
reduced market model can be defined if we consider Cauchy problems where 
a Hormandcr condition is satisfied at each state space point x with respect 
to a d- dimensional subspace where this subspace and its dimension d = d(x) 
may depend on that state space point. It is worth considering this general 
case and its relation to the case where the subspace and its dimension arc 
fixed. We shall state a generalization of our main theorem in this case. It has 
the advantage of having a coordinate-invariant formulation. Consider a matrix 
function (t,x) — > (vji) n,m (t, x), 1 < j < n, < i < m on R™, and m smooth 
vector fields 

n d 
^ = EMM)^, (41) 

where < i < m. Note that we allow for time-dependent coordinates now, but 
this is not the main point. Consider the Cauchy problem on [0,T] x R™ (where 
T > is an arbitrary finite horizon) 



at 2 

u(0,x) = f(x) 
This may be rewritten in the form 

du _ 1 * \ gp- 



(42) 



(43) 

u{0,x) = f(x), 



where 



{v* j )(t,x)='£(V i )® 2 . (44) 
fe=i 

A general reduced Cauchy problem is defined by the condition that for all t € 
[0,T] and x G R" 

(v* 3 )(t,x) (45) 

has rank d = d(t,x) < n, where for each (t,x) G [0, T] x R n the number d(t,x) 
is determined by the Hormandcr condition at (t,x) £ [0, T] x K", i.e. 

{^(i, x), [V J; Vfc] (t, x), V k ],Vi](t,x),---\l<i<m,0<j,k,l---< m} 

(46) 

spans a linear subspace Wrt,x) °f dimension d(t, x). Note that at each point 
(t,x) the linear subspace is now isomorphic to M. d but varies in R" ; i.e. the 
components of the vectors in W(t, x ) are n °t fixed components of R ra in general. 
Note that the situation is different from the situation in [9] and all the examples 
above where the family (Vi)i<i< rn spans an invariant d-dimensional subspace 
of invariant dimension as (t, x) varies. The same holds for the usual stochastic 
volatility extensions (e.g. the SABR- model listed above). Let us consider a non- 
linear transformation which diagonalizes the diffusion matrix (a*j) > locally. 
Consider a function w in transformed coordinates 

w(t,y) = u{t,x). (47) 
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Then 



^ d 2 u = a * [ d 2 'w dy k dyi | dw d 2 y k 

ij — l J ij — 1 — 1 J k J 

hence the leading order term of the PDE becomes 

Tr (J x (y) A* J x (y) T D 2 y v), (49) 

where A* is the diffusion matrix with respect to the new coordinates, i.e. 
A*(t, y) = {a*j)(t, x(y)) and J x (y) denotes the Jacobian. Let diag (A^ , • • • , Xi d ) 
denote the matrix with d diagonal elements which are not equal to zero in the 
rows 11,12, ■ ■ ■ ,id, and zero elsewhere. If we have a diagonalization of (ffiy) at 
a fixed point (t, x) with d eigenvalues strictly larger than zero (i.e. a point of 
nondcgencracy in a linear subspace of dimension d) , then the implicit function 
theorem applied to 

(t,x) -> (t,x, J x (y)A*J x (yf) = [x, diag (A n , • • ■ , X id )) (50) 

shows that there is a transformation locally around that point (t, x) where we 
get a partition of the space 

K" = {U ilii2 ,..., id \l<i l <d}UC, (51) 

where the open sets 

U iu i 2 ,~ ,i„ = {x\J x {y)A*J x (y) T = diag (X h , ■ ■ ■ , X^)} , (52) 

have Xi j ^ and C is the complentary set in K™, a set of points where elliptic 
degeneracy appears in a subspace of dimension d. Note that a uniform ellipticity 
condition can play the role of a homotopy invariant and therefore it may be 
natural to write the equation in block structure in such cases. If there are elliptic 
degeneracies as in the multivariate SABR-models considered above we may use 
the additional structure to provide the reduction to the block structure (cf. 
the discussion above). If wc consider timc-discrctization, then we may consider 
different Block structures at different time steps such that the subspace may 
move in time. This leads to quite flexible schemes which cover a considerable 
class of stochastic volatility models. However, if the reduced model is obtained 
from a full factor model by a PCA (dependent on the state space point x or 
(t, x)) we are in the general situation of a Cauchy problem (42) (and/or (43)) 
which cannot be reduced to a block structure. Let us consider the essential 
time-homogenous case (i.e. coefficients depend on the spatial variables). In this 
case we may consider the intersection of all x-dependent subspaces which are 
spanned by the local Hormander condition at x, i.e. 

Ih ■■= rixgR-Wa (53) 

where W x is the vector subspace of dimension d{x) spanned by (46) above 
at x. Here the notation Ih indicates that we are considering an intersection 
of spaces defined by local Hormander conditions. Our most general theorem 
will show that regular global solutions of (42) (rsp. (43)) exist if the data 
arc rough (i.e. in L p (R Tl ) only in Ih and are smooth on the complementary 
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vector subspace R™ \Iu- In the following section we look at the situation from 
the perspective of elementary stochastic analysis. We observe that from this 
perspective regularity is closely linked to regularity of the data. We shall see 
then that a constructive analytic scheme together with Malliavin type estimates 
lead us to stronger results. 

5 Existence of global weak solutions and regu- 
larity from the stochastic perspective 

First let us note that we use the term 'weak solution' here as it is used in the 
context of partial differential equations, i.e. equations can be interpreted in a 
weak sense and solutions can be constructed in weak function spaces. This use of 
the term 'weak solution' is quite different from its use in the context of stochastic 
differential equations. Here, a solution where the versions of the stochastic terms 
are given in advance and the constructed solution is adapted to a given filtration 
is often called a 'strong solution'. In the latter sense the solution of the cited 
theorem 3 is a strong solution but, as we observe, the solution exists only in 
a weak function space, hence is 'weak' from an analytical perspective. The 
existence of stochastic processes related to degenerate parabolic equations is 
not unfamiliar in the context of ordinary stochastic differential equations or in 
the context of more advanced stochastic analysis (such as Malliavin calculus). 
A standard theorem of ordinary stochastic differential equations (for statement 
and proof cf. [20]) is the following. 



Theorem 3. Let T > and let b : [0, T]xl" W 1 , and a : [0, T] x R™ ->■ K : 



for some constant generic C > and with \a(t,x)\ = y^2ij \ a ij\ 2 (1-1 denoting 
the Euclidean norm), and such that 



\b(t,x)-b(t,y)\ + \a(t,x)-a(t,y)\<C\x-y\; x e M", te[0,T]. (55) 



Let Z be a random variable independent of the a-algebra generated by 
W(s), s > and such that E(\Z\ 2 ) < oo. Then the stochastic differential 



dX(t) = b(t, X{t))dt + a(t, X(t))dW(t), < t < T, X(0) = Z (56) 

has a unique t-continuous solution (t,uj) X(t,uj) 7 where each component of 
X(t 7 uj) belongs to the space 

V(0,T) :={h{t,u) : [0,oo) x satisfies (i),(ii), (Hi)}, 

along with the conditions 

(i) (t, lo) — > h(t, uj) is B x J 7 -measurable, where B denotes the Borel a-algebra 
on [0, oo), 

(ii) h(t,u>) is Ft -adapted, 



be measurable functions, where 



\b(t,x)\ + \a(t,x)\ <C{t+\x\); te[0,T], 



(54) 




equation 
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(Hi) E J Q T h(t,u) 2 dt < oo. 

Note that there is no strict ellipticity condition involved. Since we are inter- 
ested in Greeks, our main concern is the regularity of expectations 

u(t,x) = E*(f(X t )), (57) 

where X t is the solution of (56) and / is some function. Regularity of u clearly 
depends on the regularity of /. What do we know from the perspective of 
stochastic analysis and what do we know from an analytic perspective? Let us 
first consider regularity from the perspective of stochastic analysis. Examination 
of the proof of theorem 3 above shows that the function (t, x) — > u(t, x) is Feller- 
continuous, i.e. is a bounded continuous function for t > if / is bounded and 
continuous. If / is non-negative and continuous, then the function u is lower 
semi-continuous, because in this case u can be obtained as a limit of an increasing 
sequence of continuous functions. However, if data are smooth then we have 
what we may call 'hypoellipticity from a stochastic perspective'. Consider the 
situation of theorem (3). If in addition the data are smooth, i.e. / <G C°°, and 
(Ty , bi £ C£° , then the solution 

(t,x)^E(f(X t )) (58) 

is smooth, i.e. in C°°. The proof is by iterative use of Dynkin's formula. 

Advanced regularity issues are discussed in the literature in the context 
of Malliavin calculus (examples of results may be found in [2] and references 
therein). Malliavin calculus in particular provides sufficient conditions for the 
existence of densities. However, as our simple example in section 1.2 shows, in 
general densities exist only in a distributional sense and the results concerning 
the existence of regular densities do not apply. Note that theorem 3 provides 
us with a solution X t with initial data x € R™ and with associated probability 
measure P such that for each t > we have a distribution 

Ff(z):=P(X t <z) (59) 

for all z = (z\, ■ ■ ■ ,z n ) £ R n . Let us consider the time-homogeneous case. A 
density (t, x; y) — > p(t, x\ y) (if it exists) satisfies 

*f(*)= / 1 ••• [ " P(t,x,y)dy (60) 



for t > and x € W 1 given. Such a 'function' p may be found in a generalized 
sense in the space of tempered distributions S' (R™). More precisely, we may 
measure it in a space 



(R n ) := {/ G S' {R n ) 1/ is a function, and ||/||2 < oo} , (61) 



11/11? = / \m\ 2 (l + M 2 Ydv, (62) 



H' 

where 



/ denotes the Fourier transform of /, and s G R. From this perspective we 
surely have 

V->p(jt,x,v)€H- n (R n ), (63) 
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but not much more. 

In analytical terms random variables are measurable functions (L°°), and the 
theorem 3 says essentially that global Lipschitz conditions of drift and volatility 
coefficient functions and initial conditions with finite second moments imply the 
existence of a global solution with finite second moments. Again, in analytical 
terms this means that the solution is a measurable L 2 -function. From an ana- 
lytical perspective, regularity of degenerate operators with regular coefficients 
(which is the 'practical perspective' from our point of view) is considered in the 
context of sufficient conditions for hypoellipticity (cf. [12]). 

6 Global regularity for a class of degenerate parabolic 
equations 

Practical incentives have led us to a class of semi-elliptic problems where the 
computation of Greeks seems to be difficult. In a situation of bounded con- 
tinuous payoffs the standard theory of stochastic differential equations provides 
Feller-continuous value functions which solve the associated Cauchy problems. 
If the payoffs are only continuous (but nonncgative) then the standard theory 
provides only lower semi-continuous solutions. For the computation of Greeks it 
is desirable to have more regular solutions. Note that in the situation of finan- 
cial applications the restriction to bounded continuous payoffs (or, equivalently, 
initial data) is a limitation. We would like to have at least Lipschitz contin- 
uous payoffs where the growth of the payoff has an exponential bound (note 
that Cauchy problems in finance are formulated in logarithmic coordinates). 
Furthermore, we have observed that there are some limitations concerning the 
regularity of the initial data. In an extreme case, where the initial data are 
(^-distributions, the solution would be a density. However, in general there is no 
density which exists in a space of i°°-functions in general. The reason is quite 
obvious: on a certain subspace the degenerate operator of the factor-reduced 
problem operates as a vector-field, and the vector-field merely transports ir- 
regularities and even distributions in a certain sense as our example in section 
1.2 shows. A way out is to have slightly different and additional assumptions 
tailored by the practical needs and by the needs of our constructive analytical 
method. We allow the initial data to be of exponential growth at infinity. On 
the other hand, they are allowed to be L p on a subspace where the diffusion 
lives. In this way we avoid transportation of the kinks by the vector field in our 
analytic AD-scheme 4 . Next we state and prove an extension of [9]. The proof 
method is similar, but in addition we use certain estimates obtained in [23]. We 
formulate the problem in a coordinate-dependent way and state a coordinate- 
independent formulation below. We state the theorem in the essential case of 
space-dependent coefficients. 

Theorem 4. Let < d < n, T > some real number, and 1 < p < oo. 
Consider the Cauchy problem (7) on [0,T] x R n . Assume that the Hormander 
condition holds on the the subspace M. d of the first d coordinates. Assume that 

4 Originally AD(I) scheme refers to an alternate direction implicit scheme for solving PDEs. 
Since the iteration scheme used in theorem 4 closely resembles this technique, we sometimes 
refer to it as an analytic AD(I) scheme. 
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(64) 



the initial data function f : R n — > R satisfies 

(i) for all Xd+i, • • • ,x„ fixed the function 

(Xl, • • • ,Xd) -> /(Xi,- • • • • • ,x n ) 

«s £f oc , 1 <_p < oo on R d , 

(m) /or all xi, ■ ■ ■ ,Xd fixed the function 

(Xd+l, ■ ■ ■ ,X n ) -> /(Xl, • • • ,3^,2^+1, • • • ,X„) 

is C°° (R n - d ), 

(Hi) for all x £ R" 

< C exp(C|x|) /or some constant C > 0. 

ylssitme </ia< i/ie coefficients are smooth (i.e. C°° ) of linear growth with bounded 
derivatives, i.e. 

a 3l g (R™) (65) 

/or i = and 1 < j < n, or 1 < i < m and 1 < j < d. Then the Cauchy 
problem (8) has a global classical solution u, where 

ueC°°((0,T]xR n ), (66) 

where for data f the singular behaviour int = 0is determined by the Malliavin- 
type estimate in [23] as follows: for given natural numbers m and N there is a 
number q such that the solution u and its time derivatives up to order m and 
its spatial derivatives up to order N are located in the space 

C?n°N ([0, T] x R n ) := {v\t.H e C m>N>loc ([0, T] x R")} , (67) 

where 

D^gWioc < oo 

[ l<m \ol\<N 

(68) 

with ||.||; oc denoting the local supremum norm. Moreover, a lower bound for q 
is given by q = max|a| < JV(n mjtl .o — n/2) where n m>a _Q is determined by the 
estimate in [23] of the singular behavior of the density (cf. theorem 13 below). 

Remark 5. Note that the space in (67) is just introduced in order to consider 
the singular behavior at t = 0. On may sharpen this estimate (considering 
subspace dimension d). 

In the proof of the theorem we construct a functional series which converges 
to the solution of the Cauchy problem. The convergence is proved by applying 
Hormandcr estimates for subproblems related to the functional scries. 

Remark 6. Note that the exponential growth of data typically appears in the 
context of finance- cf. a call option in the logarithmic coordinates. 

Next we state the corollary, which can be obtained directly from the proof 
of theorem 4. However there are some special issues like the computation of the 
time step size, where the step size may depend on the regularity of the data 
such that schemes with increasing time step size seem possible. However, these 
are special computational issues which will be considered in the second part of 
this paper. 
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Corollary 7. The solution of the Cauchy problem (8) can be obtained by com- 
putation of the functional series 

oo 

u(t, x) = U 

i=i 

where the terms of the right side of (69) are defined below. The series is com- 
puted with respect to some time discretization of the form 

pT ,2pT , . . .,mpT , . . . 

for some Tq > 0, where p can be computed a priori and depends only on the 
diffusion dimension k, the dimension n, and on the coefficient data. Assume 
that u has been computed up to time mpTo for some m > 0. Then we compute 
u on the extended time horizon [0, (m + l)pTo] as follows. We first compute the 
solution of the first order Cauchy problem 

¥-e:Ui/^)|£ = o, (to) 

with initial data u(mpTo, Xi, . . . , Xd, •) for fixed (xi, . . . , x<j). We may compute 
this solution by solving associated ODEs along characteristic curves (which are 
also solutions of ODEs). Then we compute the solution of the parabolic Cauchy 
problem 

§u)__spk a* (x) d 2 " 1 T d u (x)^- 

(71) 

= 'l2i=d+l t l i( X )'dx~ ^ 

with initial data u(mpT$, •, Xd+i, ■ ■ ■ , x n ) for fixed (xd+i, ■ ■ ■ , x n ). We may com- 
pute this solution using higher order density approximations. Alternating, we 
next compute for I > 1 the solution 

Su 21 = u 21 - u 21 - 2 (72) 

of the equation 

at Z-<i=(i+iptW dxi 



J2iJ = l a ij( x ) dxidxj + S»=lM»( a V dxi 

and the solution 

Su 2l+1 = uP- 2l+1 - u 21 - 1 (73) 

of the equation 

dSu 2l+1 * „d ,\dSu 2 



dt J2i,j = l a *j( x ) dx z dxj Z^lMifa) dx, 

V™ a (x) dSy21 

l^i=d+lt l i-\- L > dxi ' 



(74) 



where for m > 1 5u p,m has zero initial conditions, i.e. 5u p,m (0, x) = 0. Here 
for m = 1, 3, • • • (xd+i, • • • , x„) is /ixed, and /or m = 0, 2, • • • (xi, • • • , Xd) is 
fixed. 
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Proof. (Theorem 3). Note that for d < n systems of type (8) are equations 
of mixed type: one part of the equation acting on a d-dimensional subspace is 
hypoelliptic while the other part may be termed hyperbolic (a very simple form 
of scalar first order type). However, fixing x d = (xi, ■ ■ ■ , xj) we may interpret 

Y, »i(x d ,x n - d )— (75) 

i=d+l * 

as a vector field acting on the subspace R"~ d . Similarly, if we fix x n ~ d , then on 
the subspace M. d the operator in (8) acts essentially as in the parabolic equation 
on subspace R d of form 

§ " Et- =1 <j(x k , x n - k )Slt ~ E?=i * n ~ d )& = 0- (76) 

A standard AD-schcme may be used to join these heterogeneous parts of the 
whole space operator together. However, in order to get a convergent scheme we 
first need two basic assumptions: i) the existence of global solutions to certain 
vector field problems associated with the vector field (75), and ii) the existence 
of global solutions for the parabolic equations of type (76). If both assumptions 
are satisfied then we are able to present an AD-scheme which proves to be 
convergent locally with respect to time. Iteration of the scheme in time using 
the semigroup property of the operator then leads to a global scheme. Let 
us start with the first assumption, the global existence of certain vector fields. 
Fixing x d = (x\, • • • , Xd) we may interpret 

£ ^(x d ,x- d ) — (77) 

i=d+l 1 

as a vector field acting on the subspace In the scheme we have to deal 

with inhomogenous vector fields. 

The proof thus consists of three parts. The first part will consider the 
solution of the vector field for a finite time step. The second part will consider 
the solution to the parabolic subproblem. In the following sections the coupling 
of the two steps is represented by a function g below. The third part of the 
proof combines the two steps and shows that the scheme converges, constructing 
a solution on [0,T] x R™. However, before we start the argument, we should 
mention that it is sufficient to prove the theorem under the stronger assumption 
of bounded payoff /, because we can transform the problem (8) for u to a 
problem for 

u := e- d{x) u := e -V a +«M 2 u (78) 

for some a > 0, q > C 2 , and where |.| denotes the Euclidean norm. Then u 
solves an equivalent problem with identical diffusion term but transformed drift 
vector b := b — ^Vrf • aa T and an additional potential term c := c + b • Vd — 
^tr (crcr T ) D 2 d — ^ | V dcr | 2 . Here D 2 d denotes the Hessian of the function d and 
tr denotes the trace of a matrix. The following argument then holds for bounded 
payoffs, and applying the transformation above we see that it can be extended 
easily to the case where exponential growth of payoffs is allowed. 

Remark 8. Note that we may choose q > C 2 > such that the initial data 
decay exponentially as \x\ 'f oo. From now on we assume that q is chosen in this 
way. 
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6.1 Existence of the Vector Field 



First wc have 

Proposition 9. Fix x d E M. d . Assume that the conditions of theorem 1 are 
satisfied. Assume that g €E C 1 ([0,T] x R"). Then there exists a smooth global 
flow T l generated by the vector field (75) on [0,T] x W l ~ d such that the first 
order equation problem 

W = ELz+i x n - d )£- u + z"~ rf ), 

(79) 

u(0,x d ,x n - d ) = f(x d ,x n - d ), 

has the solution 

u(t,x d ,x n - d ) = f (x d ,T t x n - d ) + f g(x d , P~ s x n - d )ds. (80) 

Jo 

Proof. Note that the variables x d are fixed and serve as parameters. Consider 
the characteristic form 

n 

Xl(«,0 = £(>- ]T Mi& (81) 

i=k+l 

of the operator L = ^ - X^Ld+i/^^T' wnGr c £ = (£o, • • • , £«)■ Tlle 
surface 5* := {i = 0} has a constant normal vector (1,0, • • • ,0), hence is non- 
characteristic for the surface S, i.e. at any point z = (t, x) we have 

(1, 0, • • • , 0) i char z (L) := U ± 0|& - ^ = 1 . (82) 

Hence, basic PDE-theory tells us that the first order Cauchy problem has a 
unique local solution in a sufficiently small neighborhood of the surface S and is 
given in the form of solutions of associated ODEs along its characteristic curves. 
This leads to a solution up to a time T\. Next we may iterate the argument in 
time. Assume that this does not lead to a global solution but to a limit > 0. 
Then on the time horizon [0, Too] we have a classical solution. Moreover the 
solution has a representation on this horizon as a family of ODE-solutions along 
characteristic curves, and where the assumptions on the coefficients (65), (55) 
imply that this family of solutions is uniformly bounded up to time T x . Hence 
we may apply the first order PDE argument above again for the Cauchy problem 
with initial data St^ '•= {t = T^} and extend the solution beyond the horizon 
[0, Too]. Hence there is a unique global solution. For each x^~ d € R n_d the flow 
F t of the vector field J2i f J -i~Bir defines a characteristic curve Xq~ (t) := F t XQ~ d . 
Note that 

T l x n - d (83) 
is a solution of the homogeneous Cauchy problem 

du _ v^n ,, (rrd ™ n ~d\ d „, 

(84) 

u{0,x d ,x n - d ) = f(x d ,x n - d ), 

and then the form of the solution (80) of the inhomogenous equation follows 
from Duhamcl's principle. ■ 
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6.2 Construction of the solution via an AD-Scheme 



The notation above, which indicates that some coordinates are fixed (x d or 
x n ~ ), is a little cumbersome and so from here on we shall sometimes just write 
x instead of (x d , x n ~ d ) when it is quite clear from the context which components 
of x should be considered as fixed. Next we define a local iteration scheme 
involving global flows of the type discussed in Proposition 9 above and solutions 
of parabolic equations of the form 

1,3 — 1 J Z— 1 

with x n ~ d fixed. The natural ansatz is an AD-scheme of the following form: we 
define 

Vector Field Step: (I > 0) 

if I = (86) 



and 

Diffusion Step: (/ > 0) 



dt 

\ du 



V fc a* (x) d2u2 ' +1 T d a (x) du2 ' + 



(87) 



For each m we define u m (0,.) = /(.) and M m+1 (0, .) = /(.). Here, in equa- 
tion (87) we understand {xd+i, ■ ■ ■ ,x n ) to be fixed, and in (86) wc understand 
(xi, • • • , Xfe) to be fixed. In order to prove convergence time step by time step 
we rewrite the scheme in time-dilatation coordinates (p will be small) 

t : [0,oo) -> [0,oo), 

(88) 

t(r) = pr. 

Then we get an equation in r equivalent to (8) where the coefficients of the 
symbol of the operator become small if p is small. We have 



r/r 
and 



dt 

P, 



(89) 



u(0,x) = f(x). 

An iteration step of the scheme considered in transformed time r for some 
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time horizon [0,Tb] is then given by 

Ei=d+iPMi(a 



for / = (90) 

and 



dr ~ ' Pl2i,j=l a ij( x ) dxidxj Yui=\ PP-i{ x Y dxi 



(91) 



for Z > 0. The initial conditions are 

u p ' m (0,x) = f(x), m>0, (92) 

where for m = 1, 3, • • • (xd+i, • • • , x n ) is fixed, and for m = 0, 2, • • • (xi, • • • , x^) 
is fixed. We construct the solution in the form 



(r,x) =u"- 1 (T,x) + ^(5u' , ' 2 ' +1 (r,x), (93) 
where for Z > 1 

<Ju'' 2,+1 = uP> 2l+1 - n"' 21 - 1 (94) 



(95) 



satisfies 

dSu^l_ oT k a *.( x ) d 2 6u^+i _ T d (x) §M^l 

dr P £—fi,j=l ij V / 8xidx;,- Z^i=lPPt\ x ) dx, 

- V™ ou fxl a<5 " P ' 2 ' 
and in each substep where the right side in (95) 

5u p,2l = U PM _ u p,2l-2 (Q6) 

satisfies 

— 2^i,j=l P a ij\ x ) dxidxj 2si=\ PP'i'K ) dx t 

Moreover, for m > 1, Su p ' m has zero initial conditions, i.e. Su p ' m (0, x) = 0. 

Remark 10. The construction of the solution is designed to be used as a nu- 
merical scheme in the following sense 

1. the suitable time step size pTo may be determined explicitly, 

2. the irregularity of the payoff affects only the first summand w p ' 1 of (93). 
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We claim that for small p the scheme just described is locally convergent 
with respect to time. Then iteration of the scheme in time using the semigroup 
property leads to a convergent scheme of a global solution to the Cauchy problem 



^-|E*;=iK-(z) 
u p {0,x) = f(x). 



Er=iPMi(*)i£ = o 



(97) 



Note that iteration in time means that we start the next time step with 
the initial data u p (Tq, .), and after repeating the scheme above we get the next 
initial data u p (2Tq,.) and so on. The choice of p depends on certain a priori 
estimates which we cite below. Since the constant of the a priori estimate 
depends only on time-independent parameters we can set up a uniform time 
discretization. Hence, the solution can be computed by a time discretization 
where the original problem is computed on subintervals [ipTg, (i + l)pTo\, < 
i < M for some integer M with T = MTq and some small p. At each time 
step i this is equivalent to computing the transformed problem on subintervals 
[iT , (i + 1)T ], < i < M, where the series 

u p (t,x) =u p ' 1 {t,x)+^28u p - 21+1 {t,x), (98) 
i>i 

converges pointwise to a classical solution of (8). Since the equations are only 
hypoelliptic on the subspace R d we cannot apply Schauder estimates (at least 
not directly) . However estimates in [23] based on the Malliavin calculus lead to 
estimates of densities at each iteration step on the subspace R d . This leads to 
an estimate of the solution where classical estimates of distributional solutions 
of heat equations with source terms are involved. The form of (98) allows us to 
prove the regularity result for data which are only L p on some subspace. It is 
convenient that the higher order corrections Su m , m > 1 have zero initial data 
as we shall see below. For u 9 ' 1 we obtain interior estimates which are Holdcr 
up to the boundary using a simple trick. Consider the series 

TCI 

u p > 2m+1 {T,x) = u p ^{t,x) + s ^5u p:21+1 (t,x). (99) 

i=i 

The local convergence of the series (98) in time is such that the time derivative, 
and the spatial derivatives up to second order of (u p ' 2m+1 ) m have an absolutely 
convergent majorant. Hence we may check that (98) is a solution by applying 
the operator of (100) to (98) and evaluating it term by term. By a simple 
induction argument you check that 

The question now is in which classical space we have a convergence of the 
functional series (99) (note that we need a space such that the derivatives of the 
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functional series up to first by member). The guide is the density estimation in 
([23]). Accordingly we define for each positiv real number q 

Clf loc ([0, T] x R") := {v\t q v G C^ oc ([0, T] xl")}, (101) 

where 

C^ oc ([0,T] x K") := J /|||/||L° C + ||/ t |lL° C + E W D *ft C < 00 \ ( 102 ) 

[ |a|<2 J 

along with ||.||^ c denoting a Holder norm. In general (in order to prove higher 
regularity) we may consider the spaces 

G q m a N C ([0, T] x R») := [v\t«v G ([0, T] x R") j , (103) 

where 

<° c ([0,T]xIR"):= |/|||/||„ + X;Pt/lla+ E IP*/IU<°4- (104) 

[ i<rn |a|<iV J 

It is easy to check that for each q > the space Cf 2 "' loc ([0, T] x W 1 ) is locally 
a Banach space where the norm is obvious from its definition and denoted by 
|| . ||^2 . The considerations above imply that proof of convergence \8u p > 2m+1 \\ 2 J. 
does indeed imply that the functional series (98) is a solution of the Cauchy 
problem. Since the equation is only hypoelliptic on a subspace we cannot apply 
Schauder estimates directly. However, we have the following Malliavin estimate 
of the density which is proved in [23]. 

Theorem 11. Consider a d- dimensional diff fusion process of the form 

d d 

dx t = Ew) d *+E CT «( x *) dW ? ( 105 ) 

i=l j=l 

with X(0) = x G M. d with values in R d and on a time interval [0, T]. Assume 
that hi , (Ty G Cf£ . Then the law of the process X is absolutely continuous with 
respect to the Lebesgue measure, and the density p exists and is smooth, i.e. 

p:(0,r]xR li xl li ->Re C°° ((0, T] x R d x R d ) . (106) 

Moreover, for each nonnegative natural number j, and multiindices a, ft there 
are increasing functions of time 

Aj, a ,p,B j>a>0 : [0,T] ->R, (107) 

and functions 

nj, a ,0,m jiat p : N x N d x N d N, (108) 

such that 

Q3 g\a\ g\0\ 



dp dx a dyP 



V2 



A 1a (t)(l + x) m i^ ( „ ,Ax-y) 
P(t,x,y) < ^ - ; exp \-B j>a , (t)^-^- 

(109) 

Moreover, all functions (107) and (108) depend on the level of iteration of Lie- 
bracket iteration at which the Hormander condition becomes true. 
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We cannot apply the theorem (11) directly, but it has a probabilistic side. 
We note 

Corollary 12. In the situation of (11) above, solution Xf starting at x is in 
the standard Malliavin space D°° , and there are constants Ci^ q depending on the 
derivatives of the drift and dispersion coefficients such that for some constant 

m u <c ljq {i + \x\y^. (no) 

Here \-\i, q denotes the norm where derivatives up to order I are in L q (in the 
Malliavin sense). 

For a proof consider [23]. Using this corollary we may show that (97) and 
the right sides (95) exponentially decay as the modulus of the spatial variables 
go to infinity (recall that we consider a scheme for the transformed equation 
according to (78). We have 

Lemma 13. There is a constant c > such that for I > 1 

K'lo, \Su p '% < exp(-c|x|) as \x\ f oo. (HI) 

Here, |.|o denotes the supremum (recall that we transformed the problem along 
with its initial data.) 

Proof. First we consider the Cauchy problem for u p ^ in (97). (for each l,q) 
rewrite the probabilistic reoresentation of the Cauchy problem in terms of the 
related Markov family Y t ,x with (1 + Ixl) 7 '' 9 ^ = X t ' x , where X\' x is the 
Markov family related to the Cauchy problem (97). Then we have a probabilistic 
representation of the solution u 9 ' 1 of (97) with 

u^(t,x):=E(f(Xf)) = E(g(Yn), (112) 

and where g{.) := h((l + |a;|) 7 '".) is bounded and decays exponentially for 
\x\ t oo . Similarly for equations with a source term. The densities related to 
the Markov family Y x have densities which can be estimated by densities which 
satisfy an equation which is equivalent to the heat equation with some initial 
data or some source term. Hence there is a constant c > such that 

K^lo < exp(-c|x|) as \x\ f oo. (113) 

Next |<5u p,2 |o < exp(— c\x\) as \x\ f oo follows from proposition 9. Similarly, for 
I > 1 we use the probabilistic representation 

5uP< 2l+1 (t, x) := ft E (6uP' 2l {s, Xf)) ds 

(114) 

= /* E (SuP' 2l (s, (1 + |a:|)T'.«y/)) ds, 

and proposition 9 and get \5u p ' 2l \ , |<5m p ' 2Z+1 |o < exp(— c\x\) as \x\ f oo induc- 
tively. ■ 

Now we may write the members of (99) in terms of the family of densities of 
the Markov family (Yf)^^. First we define the family of densities (i, x, y) — > 
p Y (t,x,y) via 

P x (Yf edy)=p Y (t,x,y)dy, (115) 
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where P x denotes the law of the process Y t x which starts at x £ M™. Then we 
get the representations 

u^(t,x) :=E(f(X?)) =E(g(Y t x )) = J g(y)p Y (t,x,y)dy, (116) 

and 

SuP' 2l+1 (t, x) := /J E (5uP> 2l {s, X*)) ds 

= f* E (6uP> 2l (s, (1 + \x\)v<<Y?)) ds (117) 

= /„* 6uP> 2l (s, (1 + \y\)^"p Y (*, x, y)dyds 
From the estimates (110) (proved in in [23]) we get 

\Y?\i, q <Ci, q (118) 

for some constants < Ci. q < oo, and for the choice ji. q > too,o,o 

I ^ , A),0,o(*) ( r> / i \ (^-2/) 2 ^ , 110 x 

|py(t,x,j/)| < t » 0OO exp I -B , 0) o(*) - f ) (119) 

Since Ao.o.o(0 and i?o,o,o(i) are increasing functions of t (cf. [23]) we may 
choose A* := A^o.o^) and B* = i?o,o,o(0) in order to get the estimate 

A* f (x-v) 2 \ 

M*,*,»)| < —exp^-B'^L) (120) 

Similar estimates hold for higher time and higher spatial derivatives of the 
density py of course. The only difference is that the singular behavior with 
respect to t as t J, and the constants A* , B* change according to the estimate 
(109). More precisely, for each nonnegativc natural number j, and multiindiccs 
a, p we get constant 

A j,a,p = A j,^p{T),andB* jap = n,,,.Ji)\. 

and functions 

nj, a ,p,m jta ,p : N x N d x N d -> N, 

as in [23] such that 
Qj g\a\ g\p\ 



dP dx a dyP 



p Y (t,x,y) 



-exp \-B j<aj3 I =: q J Y ' p {t,x,y). 



(121) 

Now by Youngs inequality, for / € L and g E L p , < p < oo, we have f*g £ L p 
(* denoting convolution), where 

\\t*9\\v<\\SU\9\\p- (122) 
This can be used to estimate general functions s £ L 1 (M x K n ) 

sup T / ^(a.vJ^^a.yJIdsdy^CHalli (123) 
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for q > rij ia> — n/2 where C is an upper bound of the L\ norm of the kernel py. 
Now we may apply such estimates for source terms 6u p ' 2l+1 and their derivatives 
for I > 1 (an extension which is zero for all t < can be defined in order to have 
them in the form of s above), and similarly for u p ' . More precisely in order to 
get convergence in 



^-tq,et,loc 



([0, T] x R") := {v\t q v e C;'; c ([0, T] xl")}, (124) 



of the functional series (99) and its time derivatives up to order m and its spatial 
derivatives up to order TV we act as follows. For given Tq > and any nonnega- 
tive numbers m, N we find p small enough and q such that iterative application 
of the estimate leads to a absolutely convergent geometric series majorant for 
the solution approximations u p,m = u p:1 + Y^iLi 5u p ' 2l+1 and for their time 
derivative and their spatial derivatives up to order m and N respectively. We 
now have for some c € (0, 1) 

\^ 2l+1 \ q m a N C < c2l+ \ \Su p ' 21 \Xn C < c 2 ' I as 1 1 oo. (125) 

which implies convergence on the domain with horizon Tq for the appropriate 
choice of p. The proof of convergence for one time step uses a priori estimates 
of the Malliavin type for the equation in transformed coordinates. This corre- 
sponds to one time step on the subinterval [0,pXo] m original coordinates (here 
To > is arbitrary while p has to be chosen). It is clear then that the same argu- 
ment can be applied to all successive subintervals [ipTg, (i + l)pTo], < i < M 
resp. [iTo, (i + l)Tb], < i < M with respective initial data obtained from the 
previous time step (we use the semi-group property of the evolution equation). 
Iteration of this argument in time proves global convergence of the scheme. Note 
that for each time step the functional series and its derivatives up to first order 
with respect to time and up to second order with respect to the spatial variables 
have an absolutely convergent majorant. Hence componentwise differentiation 
of the functional series is justified, and the convergent functional series is indeed 
a solution of the Cauchy problem. ■ 

Remark 14. For numerical purposes (in order to avoid dealing with time sin- 
gularities numerically) you may define a corresponding series for the function 

(t, x) -> v p (t, x) := -t q u p (t, x). (126) 

q q 

Indeed, in order to prove that (t,x) — > v p (t,x) is j-times differ entiable with 
respect to t and a-times differ entiable with respect to x we may choose 

q ■= n>j,a,o + 1- (127) 
The iteration scheme is then of the form 



- £-=i PlH(x)S^l = E t k+1 PIH{x)9%£- + pt^Su^it, x), 
with zero initial conditions, i.e. Sv p,2l+1 (0,x) = 0. 



(128) 
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7 A coordinate independent generalization of Hormander's 
result 



We state our most general theorem in the essential time-homogenous case. Con- 
sider a matrix- valued function x — > (vji) n - m (x), 1 < j < n, < i < m on M. n , 
and m smooth vector fields 

n ^ 

^^(x) — , (129) 
3=1 3 

where < i < m. Note that we now allow for time-dependent coordinates, but 
this is not the main point. Consider the Cauchy problem on [0,T] x M. n (where 
T > is an arbitrary finite horizon) 

(130) 

u(0,x) = f{x). 
Define for all x E R n 

W x :=span{ ^(x), [V,-, 14] (x), 

(131) 

[[Vj,V k ] , V{\ (z), ■ • • |1 < i < m, < i, k, I ■ ■ ■ < m}. 

Furthermore let Ih = n x ^M"W x - The following general theorem assumes that 
the roughness of the data is located in Ih- Note that in the case Ih = {0} the 
following theorem collapses to proposition (??) above. We have 

Theorem 15. Let 1 < p < oo. Consider the Cauchy problem (130) on [0,T] x 
R™. Assume that the initial data function f : W 1 — > K satisfies 

(i) the function f is £f oc , 1 < p < oo on Ih, 

(ii) the function f is C°° on R" \ Ih, 

(Hi) for all x E K" 

1/(^)1 ^ Cexp(C|x|) for some constant C > 0. 



(132) 



Assume that the coefficients are smooth (i.e. C°° ) and of linear growth with 
hounded derivatives, i.e. 

va E C% (R n ) (133) 

for i = and 1 < j < n, or 1 < i < m and 1 < j < n. Then the Cauchy 
problem (8) has a global classical solution u, where 

u E C°° ((0, T] x W 1 ) , (134) 

where for sufficiently regular data f the singular behaviour ist = 0is determined 
by the Malliavin-type estimate in [23] as follows: for given natural numbers m 
and N there is a number q such that the solution u and its time derivatives up 
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to order m and its spatial derivatives up to order N are located in the space ( we 
do not consider a Banach space here, cf. remark to 4 above) 

C q m a i° C ([0, T] x K") := [v\t«v E C%% c ([0, T] x R")| . (135) 

Moreover, q = max|a| < N (n m , a ,o — n/2) where n miQj o is determined by the 
estimate in [23] of the singular behavior of the density (cf. theorem 13 below). 

The proof for this theorem is a quite similar to the proof in the preceding 
section. The difference is that instead of the flow of a vector field we now have to 
consider a flow (defined by the semigroup) of a semi-elliptic SDE in a subspace 
of smooth data. 

8 Higher order probabilistic weighted Monte Carlo 
schemes and comparison with proxy schemes 
and partial proxy schemes 

The abstract scheme above does not show us how to compute. It especially does 
not show us how to compute in a probabilistic manner, i.e. with Monte-Carlo 
schemes. In this section we present weak higher order weighted Montc-Carlo- 
schemes in a way that allows us to extend the scheme of [10] and related schemes 
in [7] and [18]. Why do weak higher order schemes matter in finance (compared 
to Euler-schemes, Milstein schemes and other classical schemes)? 

In order to answer the question let us recall some basic facts. Consider an 
ordinary stochastic differential equation of the form 

dX t = b{t, X t )dt + a{t, X t )dW t , X = x e R" (136) 

on a time interval [0, T] and a time discretization fcAT, < k < N with 
NAT = T. We assume that the assumptions of theorem 3 above are satisfied. 
An Euler scheme is described by the recursion 

Y k+X = Y k + b(kAt, Y k )AT + a(kAT, Y k ) (W {k+1)AT - W kAT ) . (137) 

The latter scheme is a (simple) example of strongly convergent scheme, there is 
some 7 > such that 

e(<5) = E(\Xt — Y(T)\) < CAT 1 (138) 

as AT becomes small. Indeed, in the case of an Euler scheme we have con- 
vergence with 7 = 0.5. More elaborate schemes like the Milstein scheme and 
trapezoidal methods also converge in a strong sense. Why should we consider 
higher order weak schemes then? What is the relation between strong conver- 
gence and weak convergence? Weak convergence is defined with respect to a 
function class. We say that a an approximation scheme Y converges weakly 
with order j > to X as At 4- an d with respect to a function class C, if for 
all g G C 

\E (g(X T )) - E (g(Y T )) \ < CAT 1 (139) 

as AT | 0. It depends on the regularity of functions in C whether (139) is a 
strong condition. However, especially in finance, we may have low regularity of 
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payoffs, e.g. in the case of digital payoffs. If we compute As near maturity then 
the situation is even worse, i.e. we get distributional behavior at maturity. In 
particular, the asymptotic behavior of the density of X is not correctly approx- 
imated by an Euler scheme, because the Euler scheme leads to an Euclidean 
limit while 

lim o ATmp(AT,2-, y) = d 2 (x,y), (140) 

where with {aa T ) 1 ^ being the entries of the inverse of <jo t , the function d 2 
is the Riemannian metric induced by the line element ^™ j =1 (ccr T ) l: ' dxidxj . 
However, (140) is also interesting because it is the leading order term of analytic 
higher order expansions of the density. Numerically stable higher order density 
approximations of the form 

1 ( d 2 R (x lV y 



p(f,x;0,y) = -^exp(^S^) ( £ d k (x, y)t k ] (141) 



\fc=o / 
have been developed in [17] recently. 
Remark 16. Note the difference to the W r KB '-expansion 

p(t,x;0,y) = -^e^(-^^+f^c k (x,y)A , (142) 

used in [18]; in the latter article the expansion was applied to classical LIBOR 
market models and is computed up to c±. 

In [18] standard interest rate options with a maturity of ten years could 
be computed in the standard full-factor Libor market model framework in just 
one time step. However, there are disadvantages to simulating the full factor 
model. First, computational experience shows that more than 7 or 8 Brown- 
ian motion factors lead to numerical noise and do not improve the calibration 
of models with respect to the market data. Furthermore, there is a trade-off 
between the efficiency gained by the analytic higher order expansions and the 
loss of efficiency caused by the number of Brownian motion factors of the full 
factor model. Naturally questions arise. Is it possible to extend the use of 
higher order analytic expansions in the framework of the proxy scheme in [10] 
and/or in the partial proxy scheme in [7]? Second, can the latter schemes be ex- 
tended by a probabilistic scheme defined along the constructive solution scheme 
of semi-elliptic diffusions defined in this paper. Clearly, the regularity result 
obtained in this paper is an advantage: it allows us to get error estimates for 
weak higher order schemes for Greeks which cannot be obtained without such a 
result. Furthermore, it allows us to derive estimators which are of bounded vari- 
ance even for small time-to-maturity and/or small volatilities. In this respect 
the result of this paper allows the extension of the result of bounded variance 
in [18]. Let us go deeper, and derive a probabilistic scheme from the series of 
parabolic Cauchy problems which define the solution of the semi-elliptic equa- 
tion of theorem 4 above. For its implementation we have to select a cutoff and 
approximate on a finite domain il C R™, but this is not essential for the follow- 
ing discussion. Furthermore we consider only the initial time step. It is clear 
then how to proceed in time with respect to a time discretization. For each 
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time step we assume that the diffusion coefficients a,ij and the drift coefficients 
depend oniy on the spatial variables. Let k < d be fixed and fix the coordinates 
x n = (xk+i, • • • , x n ). Let p x n-k be the fundamental solution of 

du p 



dr ' PJ2i,j=l a ij( X )dx z dx 3 

(143) 

Then according to our results above we get the probabilistic representation 
U p>\t, .) = J Rk f(y k , x n - k ) Pxn - k (T, x k ;0, y k )dy k 

+ Jo /r* EIU+i Pi(t(s),y k , x n ~ k ) x (144) 
§£-(y k , F T x n - k )p x n- k (T, x k ; s, y k )dy k ds. 

Similarly for the higher order terms we get the recursive probabilistic represen- 
tation 

Su^+Ht, .) = /; J Rk Er=, + i P ^(y k ,x- k )x 

(£- Jo Et'=i pa^x^q^iy^^x^ds 

(145) 

+£- Jo TLi pp*(y k > ^- fc )^£ zi ( s > y k > F T - s * n - k )) x 

p x n-k(r, x k ;s,y k )dy k ds. 
The representation (145) together with the initial representation (144) gives the 

k 

following 'naive' weighted Monte-Carlo scheme. Let ( x be a random variable 
with normal density <fi(t, x k , -). 5 Then we have the probabilistic representations 



U P,1 (T )_ E I f(C ^-K„- fc (r^;0,r ) 

+E^f T ZL k+ iPP>(C\x n - k )x (146) 

df_(fx k VT rr n-k\ P^-k{_r,x k -s,C k ) \ 



5 Note that in general practice <j> will be time-homogenous. 
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and 

(147) 

p..- fc (r,a*ia,C"*) dj Y 

Substitution of expectation £ by the sum -^j Em=i an< ^ writing for each m 
the m-th realization of as mC 31 leads to the following 'naive' Monte-Carlo 
scheme for the price 

V T i ■) - M Z^m=l I 0(T,x", m C aJ ' S ) j 



1/ Em=l ( Jo Si=fc+1 P^i{mC X iX n 



(148) 



a Xj lS, m l, ,J- X ) ^ S)X fc jmf a,«>) 



and 



(149) 



Note that we have the same density for both the lower and the higher order 
terms. The estimators for the Greeks are easily obtained by differentiating 
(146) and (145) and then writing the MC approximations for the corresponding 
differentiated equations. We call the estimate for derivatives 'naive' because 
it may have unbounded variance for small time and/or small variance. We 
consider two variations of bounded variance estimators. In one case there is a 
proof of bounded variance in the case of full factor models (cf. ([18])), and some 
numerical evidence which supports the theorem. In the other case there is some 
experimental evidence (cf. [7] for a different estimator which is simpler and 
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works in many cases). In this section we shall extend the proof in the former 
case. The idea of the proof is as follows. We consider the construction of the 
solution for the Cauchy problem in the form 

u p (t,x) =u p - 1 (T,x) + ^ SuP ' 2l+1 (r,x), (150) 
l>i 

where « p,1 (t, x) is defined as in (144) and Su p ' 21+1 (t, x) is defined as in (145) 
(the probabilistic expressions of the last section) . We observed that convergence 
is such that derivatives of the series can be taken term by term. This means 
that we can estimate the variance of the Greeks term by term, i.e. up to the 
regularity constraints for the coefficients we have for any multiindex f3 

For simplicity let us consider the As'. The 'naive' estimator is constructed from 
JLi,pA( t ) - F ( 9 f(C k ^~ k )p^- k (T,x k ;0,C k ) \ 

Cte," V> V — I dXl ^(T^fc.C^) / 



-£Ur£££=*+iPMC*V-*)x (152) 



df l/-x k T T„n-k\ V^-k{T,x k ;s,C ) ,\ 

^ x > 4>(s,x k x* k ) r 

(note that the derivative of / is with respect to Xi with i > k + 1, where it exists 
in a classical sense according to our assumptions) and from 

s -5W^{r, x) = e(&- C £" x PIMiC" , - n " fc ) x 



c).vi 



(153) 



(£- Jo Et =1 P^ic" , ^ . r- s * n - k )ds 

+£- Jo Eti PPiiC" , x n - k )^§^(s, C" , J^-'x^")) x 

(f)(s : X k ,£ x " ) / 

Now we can adapt the estimator for the full factor model (cf. [18]) to the reduced 
factor case. Note that we continue to work in logarithmic coordinates (for 
financial applications the following may be rewritten in lognormal coordinates 
if this is desirable). The idea is the following. In order to get an estimator 
of bounded variance even for small volatility and/or small time we consider a 
smooth function with \dg(t, x, z)/dz\ ^ 0, and an 

Revalued random variable £ on some probability space with density At where 
Xt(z) 7^ for all z and t (for example, the standard normal density). Then 
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we write £ x := g(r,x k ,^) where we assume that it has a density <f>{r, x k ,-). 
Hence the first order term of our constructive scheme of the estimator becomes 
(analytic form) 

_d_ u P,l( T \=E\ 9 f(9( T ,x k ,t),x™- k )p :crl - k (T,x k ;0,g(T,x k ,O) \ 
dxi V ' v I Q Xl cp(r,x k ,g(r,x k ,£)) I 

+e J r EIU+i pMt, x k ,z),x n ~ k )x (154) 

df_(„ „(_ „fe t\ ■pT„n-k\ P^-k(T,x k ;s,g(t,x k ,£)) , \ 
da H y b >y\ T > J: 'Sj'-r x I <f,{s,x k ,g(r,x k £)) MA I " 

It turns out that a well chosen g in (154) is sufficient to get bounded variance 
estimators, i.e. we could keep the structure of the 'naive estimator' for the 
higher order terms. For the higer horder terns we may use the naive weigthed 
MC algorithm (without contro function) since the data (of the source terms) are 
smooth and have bounded derivatives (hence we may use (153) in this case). 
Again, substitution of expectation E by the sum Y2 m =i an< ^ ^ or eacn m 
writing the m-th realization of £ as m £ gives the corresponding Monte Carlo 
estimator 

d u P,l(-r -) _ I P ^3{r,x k 



Ht ) = — V 



dxi <j>(T,x k ,g(t,x k , m £)) 



jjEfLi ltJzY,U + iPMs,x\ m O,x^)x (155) 



df („ n („ T k A x-r n-k\ P^-k{r,x k ;s,g(s,x k , m g)) \ 
Ssci V 6 >iA 6 ' x imSj:^ * ) tt>(s,x k ,g(s,x k , m £)) I' 

and (149) for the higher order corrections (in prcatice it is sufficient to compute 
the latter sum for the lowest correction term 1 = 1). The choice of the control 
function g is discussed in [18]. It can be proved that the latter estimator is 
of bounded variance even when As are to be computed in the framework of 
stochastic volatility models and volatility and/or time to maturity is small. 
It is clear that the scheme above defines an extension of the scheme in [10]. 
Furthermore, extension of the scheme to general reduced diffusion models is 
possible. The scheme can be combined with the partial proxy scheme in [7]. 



9 Further remarks 

Though there exists no regular density for the class of semi-elliptic equations 
considered in this paper, the analytical scheme proposed leads to a probabilistic 
scheme involving regular densities living on subspaces at each state of the con- 
struction. This leads naturally to weighted Monte-Carlo schemes which improve 
the accuracy of the schemes in [7, 10] . Compared with the scheme considered in 
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[18], the advantage of such schemes is that we do not need to simulate the full- 
factor model. However, the results for bounded variance of the estimators in [18] 
can be extended to the scheme above. Note that in [18] regularity of the density 
was used in order to obtain a class of estimators for low volatilities and time 
beyond the class considered in [5]. In [18] it was shown that analytical expan- 
sions of the density may be useful in order to improve the efficiency of a scheme. 
Such density expansions may also be used for the present scheme. Similar re- 
marks apply to the use of expansions of characteristic functions as in [16], of 
course. Note that in the context of stochastic volatility models sometimes only 
the characteristic function is known. This is sometimes due to additional de- 
generacies of the diffusion coefficients. The extension of our analytical schemes 
to such models and models with jump measures is of interest. It is clear that 
such an extension hinges on a priori estimates of related degenerate diffusion 
and partial integral-differential equations, and is possible in some cases at least. 
It may also be of interest to combine the present scheme with a front fixing 
scheme for American derivatives where the obstacle problem is transformed to 
a hypcrcube (cf. [15]). Numerical analysis and simulations in the context of the 
LIBOR market model will be considered in part II of this paper. 
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